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Abstract. - We study the ground state ordering and interactions between two two-dimensional 
Wigner crystals on neutralizing charged plates by means of computer simulation. We consider 
crystals formed by (i) point-like charges and (ii) charged dimers, which mimic the screening of 
charged surfaces by elongated multivalent ions such as aspherical globular proteins, charged den- 
drimers or short stiff polyelectrolytes. Both systems, with point-like and dimeric ions, display five 
distinct crystalline phases on increasing the interlayer distance. In addition to alteration of trans- 
lational ordering within the bilayer, the phase transitions in the dimeric system are characterized 
by alteration of orientational ordering of the ions. 
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The interaction between electric double layers attracted 
much attention in the past twenty years in particular 
due to the rediscovered role of ionic correlations. It has 
been known for long time that two similarly and strongly 
charged plates can attract each other in the presence of 
multivalent counterions. This has been seen in Monte 
Carlo simulations [1-3] , observed experimentally with the 
surface force apparatus [4], deduced from the scattering 
experiments on laponite dispersions [5], phase diagrams 
of charged lamellar systems [6,7], clay platelets, etc. (see 
[8-10] for a review). 

As the attraction appears on increasing counterion cor- 
relations, the phenomenon can be conveniently character- 
ized by an electrostatic coupling parameter 5 = 2q 3 l^a s , 
which depends on the Bjerrum length lg, the counterion 
valency q, and the surface charge density a s . Several 
successful descriptions of the attraction have been built 
for the strong coupling limit E — ► oo, where the correla- 
tions are so strong that the ions form a Wigner crystal 
(WC) [11] or at least a strongly correlated liquid (SCL) at 
the charged colloid surface [12, 13]. 

Although the phase diagram of a bilayer Wigner crys- 
tal has been known for some time both at the ground 
state and at finite temperature [14-18], the interaction 
between the crystalline double layers as a function of 
their separation has not been discussed in detail. The 
mono- and bilayer Wigner crystal structures have been ad- 
dressed in literature in relation to electrons above the sur- 
face of liquid helium, two-dimensional semiconductor het- 
erostructures, Mott insulators, dusty plasmas, laser-beam- 
cooled trapped-ion plasmas, or dislocation-mediated melt- 
ing transitions [19,20]. In contrast, in the soft-matter 
and biological literature the crystal structure is typically 



regarded as an auxiliary question for evaluating the in- 
teractions between the surfaces that host these Wigner 
crystals [21-23]. Typically, only the staggered hexagonal 
crystal, which wins at large distances, is considered. A few 
recent publications discuss interaction effects such as plas- 
mon oscillations using perturbation schemes with respect 
to the ground state of the double layers at large distances 
(the staggered hexagonal phase) [21,22]. 

A characteristic feature that differentiates Wigner crys- 
tals in soft matter systems from the low-temperature elec- 
tronic ones stems from the nature of the ions: The ions 
have to be multivalent to be well ordered in aqueous 
dispersions. The high ion valency also implies that the 
Coulomb contribution dominates the free energy and in- 
layer thermal fluctuations are negligible. A situation close 
to the low-temperature behaviour can be obtained, for ex- 
ample, when polyelectrolyte molecules, globular proteins 
or charged dendrimers play the role of counterions [24-26] . 
In the general case, such ions are neither point-like nor 
spherical and their shape might influence the interac- 
tion between the bilayers. For example, in the limit of 
very long ions (DNA chains between lipid bilayers) their 
anisotropy leads to orientational ordering and formation 
of two-dimensional smectic phases [27,28]. In this work, 
we present an accurate numerical solution of the ground 
state problem both for point-like (spherical) and small but 
elongated (non-spherical) counterions. We first address 
the interaction between the bilayers with point-like ions 
and then use it as a reference system to study the effect of 
ion size and charge distribution on the bilayer interaction. 

We envision a situation where the interaction energy for 
two flat parallel Wigner crystals on a neutralizing back- 
ground is measured as a function of the interlayer distance. 
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Fig. 1: Top and side view of three 2D configurations of coun- 
terions on two parallel charged plates. The open and filled 
symbols designate the ions located on the opposite surfaces. 
Set I corresponds to a single hexagonal lattice, Set III to a 
superposition of two square lattices, and Set V corresponds to 
two staggered hexagonal lattices (the notation is explained in 
the text). 



For point-like ions the scenario includes two obvious lim- 
iting cases: at large separations, the crystals do not feel 
any discreteness of each other so that each of them has a 
2D plane-filling hexagonal symmetry. At the smallest sep- 
aration between them, where both of them lie in the same 
plane, a single hexagonal crystal is formed. At a finer res- 
olution, the transition between these two phases gives rise 
to the following sequence of structures on increasing the 
interlayer distance: a monolayer hexagonal lattice (I), a 
staggered rectangular lattice (II), a staggered square lat- 
tice (III), a staggered rhombic lattice (IV), and a staggered 
hexagonal lattice (V) [14]. Here, we evaluate the ground 
state interaction by minimizing the energy at each gap 
width and in addition, consider three particular phases in 
more detail. Namely, we consider three setups, for which 
the positions of ions on one plane fall in the geometrical 
centres of the primitive cells of the identical lattice on the 
other plane: the hexagonal monolayer structure (I); stag- 
gered square lattice (III); and staggered hexagonal lattice 
(V). These three configurations are shown in Fig. 1. They 
represent the three rigid lattices on the bilayer phase dia- 
gram, meaning that their structure does not change within 
their region of stability. The intermediate structures (II 
and IV) are " soft" , so that the aspect ratio of their primary 
cell is varying with the interlayer distance [14]. Finally, to 
study the effect of the ion shape we replace each ion by a 
dimer consisting of two identical charges connected by a 
stiff spring. The dimer is allowed to rotate and translate 
in plane and two different spring lengths are studied. 

The energy of a Wigner crystal on a neutralizing 
charged plane can be written as 
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where a s is the surface charge density of the plate, q the 
counterion charge, N the number of counterions, e the di- 
electric permittivity of the medium, and R the ion-plate 
contact distance, or the ion radius for hard sphere ions. 



The characteristic lengthscale of the Wigner crystal on 
a neutralizing background is the mean lateral distance 
between the ions a±, which is defined by tt(ci±/2) 2 = 
q/a s . The energies can be conveniently expressed in 
terms of the average ion-ion Coulomb energy q 2 /(ea±): 
U = Uea±/q 2 = 4Ue/(Tra±qa s ). In this rescaled form, the 
result would explicitly depend on neither ion valency nor 
the surface charge density. If we also suppose that the 
ions form a perfect crystal so that each ion has exactly 
the same environment, one summation over all ions can 
be performed right away. Thus, we get for the energy per 
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For two such plates, in case they are parallel to each 
other and have identical ion configurations, the total en- 
ergy will contain the self-energy, 2Lq, plus the interac- 
tion terms: ions opposite plate and plate-plate. The 
ion-plate and plate-plate energies in the case of charge 
neutrality will compensate each other as the sum of dis- 
tances from each ion of charge q to the two plates is always 
h + 2R, which is exactly equal to the distance between the 
corresponding surface element of the same net charge and 
the opposite plate. Then the only interesting contribution 
comes from the summation over the ions 
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where h the distance between the ion-containing planes, 
thus giving the distance between the plates h + 2R, and 
a the displacement vector of the lattice on one plate with 
respect to the other one. The indices j and k now mean 
a summation over the ions on different plates. As the 
interlayer interaction leads to structural transformation 
within the layer, these two terms remain interconnected 
and should always be considered simultaneously. For infi- 
nite plates, the total Coulomb energy has to be calculated 
numerically. We modeled a piece of 2D crystal of size 
40 x 40 to 80 x 80 ions with lateral periodic boundary 
conditions, and calculated the energies with relative un- 
certainty of 10~ 6 . The plates were supposed to be homo- 
geneously charged. An MMM2D summation scheme [29] 
and an MD simulation package ESPResSo v. 1.9 were 
used [30]. 

We first consider the relative energies of the preformed 
lattices. Figure 2 shows the potential energy of interaction 
between the two layers as a function of the gap width h for 
the three setups. The total energy, U = 2J7 1 + {/ 12 , in all of 
them is strongly negative and decreases in absolute value 
with the gap width. The limiting value at h —* oo on each 
curve corresponds to twice the energy of a single Wigner 
crystal U\ . To calculate the true ground state energy we 
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Fig. 2: Rescaled potential energy per ion of two parallel Wigner 
crystals on oppositely charged plates as a function of normal- 
ized distance h/a± between the crystals for the three different 
counterion arrangements. The optimized structure is obtained 
by the energy minimization using a Brownian dynamics simu- 
lation at zero temperature. The transition points as reported 
in Ref. [14] are marked by filled circles. The regions of stabil- 
ity of various crystalline structures are delimited by the circles 
and marked next to the optimized energy curve. 

perform molecular dynamics simulation of two ionic lay- 
ers with the ions constrained to their corresponding planes 
but allowed to move within the plane. For each interlayer 
distance we start simulation from each of the three pre- 
formed lattices and then record the minimal energy. We 
see that the minimal energy from the free configuration 
coincides with energy of the favorable arrangement in the 
appropriate range of distances and gets lower in two in- 
termediate regions close to the transition points I — > III 
and III — ► V (Fig. 2). In these regions, one should ex- 
pect an appearance of the staggered rectangular lattice 
and the staggered rhombic one, respectively. At h = we 
obtain the rescaled energy per ion f7 12 = —3.47502, which 
is close to value —3.47493 reported in Ref. [14] (the scal- 
ing factor in our work differs from Ref. [14] by y/n). The 
single hexagonal lattice is stable in a very small region 
close to h/a± = 0. The location of the further transi- 
tion points can be estimated from comparing the energies 
of the corresponding phases. The intersection point of 
the curves I and III in Fig. 2 corresponds to a transi- 
tion into staggered square phase. The energy of the opti- 
mized structure becomes lower than that of the phase III 
at h/a± ~ 0.55, which indicates a transition into phase 
IV. The structure V prevails at h/a± > 0.77. We note 
that the transition in our work is observed at a higher 
h/a± than it was reported in [14, 18]. The interaction en- 
ergy of the two layers is shown in Fig. 3a. The energy 
curve is smooth between the transition points. The force 
curve shows a jump at h/a± » 0.77, which is expected 
due to discontinuous character of the transition between 



phases IV and V [14]. The initial decay of the interaction 
energy is close to linear, while the asymptotic behaviour 
at h/aj_ > 1 is clearly exponential. The observed be- 
haviour is close to the asymptotic decay expected at large 
distances Ui 2 {h) oc exp(— h/(a±/2w)) [12,14] (see Fig. 2). 

We now look at the interaction between crystals formed 
by Coulombic dimers for two different lateral ion sizes L = 
0.5R — 0.l7a± and L = R = 0.34a^. As an additional 
degree of freedom is involved, we expect a richer phase 
behaviour in this case. In the limit L = the behaviour 
of the system with point-like ions is recovered. Figure 
3b presents the interaction energies as a function of the 
gap width. Simulation snapshots for three different values 
of h/a± = 0.1,0.33,0.5,0.75 and 1.2, which correspond 
to different crystalline structures, are presented in Figure 
4. We see that at small and large relative distances the 
dimers in the two layers tend to be aligned while at the 
intermediate separations the dimers in each layer tend to 
orient perpendicular to their nearest neighbors in the 2D 
projection of the lattice (either same plane neighbors or 
the neighbors in the opposite plane). 

A variation of the dimer length affects only the ion-ion 
part of the total energy of the bilayer. A system with 
longer dimers has a higher self-energy of each layer while 
the corresponding interlayer part depresses the interac- 
tion. We see from Fig. 3b that the interaction between 
the layers indeed becomes weaker on increasing the dimer 
length. Roughly, the characteristic length a± decreases by 
L/2 as compared to the point charges. Then, the interac- 
tion energy would decay as U(h) oc exp(2nh/ (a± + L/2)). 
A fit to the calculated interaction energies reveals the de- 
cay length change from 0.19aj_ for the point-like ions to 
0.15a^ for the system with L = 1 (or 0.34aj_). Ultimately, 
at L — aj_ at small separations one should recover the bi- 
layer with a'j_ = a±/\/2. 

In simulations with short dimers with L = 0.17a^ at 
various gap thicknesses we observe: (i) < h/a± < 0.27 
a staggered parallelogrammetic lattice, parallel orienta- 
tion of neighboring dimers, soft; (ii) 0.27 < h/a± < 0.8 
staggered square lattice, perpendicular dimer orientation, 
rigid; (hi) 0.8 < h/a± staggered rhombic lattice, parallel 
dimer orientation, rigid. At L = 0.34q,_l, we have the fol- 
lowing sequence of structures: (i) < h/a± < 0.17 a stag- 
gered rectangular lattice, parallel orientation of neighbor- 
ing dimers, rigid, the longer lattice constant ai is roughly 
ai + L (Fig. 4a); (ii) 0.17 < h/aj_ < 0.40 we observe 
domains of rhombic lattice, rotated with respect to the 
neighbouring domains (Fig. 4b); (hi) 0.40 < h/a± < 0.67 
staggered square lattice, perpendicular dimer orientation, 
rigid (Fig. 4c); (iv) 0.67 < h/a± < 0.91 staggered parallel- 
ogrammetic lattice, parallel dimer orientation, rigid (Fig. 
4d); (v) 0.91 < h/a± staggered triangular (a honeycomb- 
like) lattice, parallel, rigid (Fig. 4e). We note that in 
contrast with some phases of point ions (the " soft" rectan- 
gular (II) and rhombic (IV) phases), the observed lattices 
in the system with long dimers (L = 0.34a^) are rigid, i.e. 
the aspect ratio of the primitive cell stays constant within 
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Fig. 3: a) Absolute value of rescaled attraction energy and 
force between two charged plates with crystalline arrangement 
of point-like ions. The dashed line shows the predicted asymp- 
totic behaviour at large h/a± [14]. b) The rescaled interac- 
tion energy between two parallel Wigner crystals formed by 
charged dimers of the indicated length in comparison with that 
for point-like ions. The dimer charge as well as the number 
density are matching those for the system with point-like ions. 



the region of stability of the corresponding lattice. 

The most interesting observation for the dimeric sys- 
tems is related to the ability of dimers to adjust the orien- 
tation to minimize the electrostatic energy. The effect is 
strongest for the square lattices where we find perpendic- 
ular orientation of the neighboring dimers. A similar phe- 
nomenon of altering orientational ordering was discussed 
recently for a one-dimensional stack of rod-like ions [31]. 
In contrast to our observation for 2D layers, the ground 
state in a staggered ID system is represented by perpen- 
dicular orientation of the neighboring ions, which then 
changes to a twisted chiral phase on increasing density. In 
our system, the reorientation of the ions happens suddenly 
on varying the separation distance between the surfaces, 
which might find an application in nano-structuring. This 
structural transformation can be best characterized by the 
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Fig. 4: Snapshots from simulations of dimeric bilayers with 
dimer length L = 0.34aj_ at different gap widths: h/a± = 
0.1 (a), 0.33(b), 0.5(c), 0.75(d), 1.2(e). The different colors 
correspond to ions located at different planes. 



scalar order parameter S = | (3 cos 2 9—1) (the second mo- 
ment of the orientational distribution function), where 9 
is the angle between the molecule orientation and the di- 
rector, and average over all dimers is taken. This value is 
plotted in Fig. 5. One can see that the phases at short and 
long interlayer distances have S = 1, which corresponds 
to ideally aligned dimers (Fig. 4a,d,e). Further on, the S 
values of 0.25 correspond to a coexistence of two preferred 
dimer orientations (Fig. 4c), which are perpendicular to 
each other. The arrangements with two preferred orienta- 
tions are observed in the region of stability of the square 
lattice. From the plot shown in Fig. 5, we see that the 
onset of the perpendicular dimer orientation as well as the 
return into the aligned state happens in the system with 
longer dimers at smaller distances. 

Finally, we note that the presented sequence of orien- 
tational transitions exists also at finite temperatures. A 
simulation performed for L/a± = 0.34 at S = 50000 and 
3 = 10000 still showed three regions with high, then low, 
and again high order parameter on increasing the inter- 
layer distance, although the short-distance and the long- 
distance phases become less aligned due to thermal fluc- 
tuations. The observed order parameter was S ~ 0.9 at 
small interlayer separations, and S sa 0.8 (S = 50000) 
and S ~ 0.4 (S = 10000) at large separations (Fig. 5, 
triangles), while the values for the perpendicular dimer 
orientations, S « 0.25, did not change with temperature. 
One can expect that the temperature region of stability 
of the square lattice with the perpendicular dimer orien- 
tation is larger than that for the aligned phases, which 
follows from the increasing stability of the square lattice 
itself in the system with point ions [16]. We also note 
that the bilayer system is stable with respect to normal 
fluctuations, which are suppressed in our setup by the re- 
pulsion between the mobile ions that confines the ions to 
thin layers at the walls. 

We calculated the ground state interaction energy of two 
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Fig. 5: Orientational order parameter for dimeric Wigner bilay- 
ers with two different dimer sizes. The value S = 1 corresponds 
to perfectly aligned dimers, S — 0.25 to a coexistence of two 
perpendicular dimer orientations. The lines are guide to the 
eye. 



planar ionic double layers represented by Wigner crystals 
of point-like charges or charged dimers, where the latter 
model mimics the situation of screening the interaction 
between charged surfaces with polyelectrolytes or non- 
spherical molecules. The crystalline structure observed 
with point-like charges agrees with that reported in the 
earlier literature, while novel structures with unusual ori- 
entational ordering were observed in the system with elon- 
gated ions. In all cases we found exponential asymptotic 
decay of the correlation attraction at interlayer distances 
larger than the characteristic lateral separation between 
the counterions. We found that the elongated ions pro- 
duce in general weaker correlation attraction attraction 
than the point-like or spherical ions of the same charge. 
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